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Abstract 

It is our intention to provide via fractional calculus a generalization of 
the pure and compound Poisson processes, which are known to play a 
fundamental role in renewal theory, without and with reward, respectively. 
We first recall the basic renewal theory including its fundamental concepts 
like waiting time between events, the survival probability, the counting 
function. If the waiting time is exponentially distributed we have a Poisson 
process, which is Markovian. However, other waiting time distributions are 
also relevant in applications, in particular such ones with a fat tail caused by 
a power law decay of its density. In this context we analyze a non-Markovian 
renewal process with a waiting time distribution described by the Mittag- 
Leffler function. This distribution, containing the exponential as particular 
case, is shown to play a fundamental role in the infinite thinning procedure 
of a generic renewal process governed by a power-asymptotic waiting time. 
We then consider the renewal theory with reward that implies a random 
walk subordinated to a renewal process. 

MSC 2000: 26A33, 33E12, 44A10, 44A35, 60G50, 60G55, 60J05, 60K05. 
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1 Essentials of renewal theory 



The concept of renewal process has been developed as a stochastic model 
for describing the class of counting processes for which the times between 
successive events are independent identically distributed (iid) non-negative 
random variables, obeying a given probability law. These times are referred 
to as waiting times or inter-arrival times. For more details see e.g. the 
classical treatises by Khintchine [12| . Cox [2J, Gnedenko & Kovalenko [BJ, 
Feller [5), and the recent book by Ross [19]. 

For a renewal process having waiting times T\ , T2, . . . , let 

k 

to = 0, t k = Y / T j , k>l. (1.1) 

3=1 

That is t\ = T\ is the time of the first renewal, ti = T% + T2 is the time of 
the second renewal and so on. In general t k denotes the A;th renewal. 

The process is specified if we know the probability law for the waiting 
times. In this respect we introduce the probability density function (pdf) 
4>(t) and the (cumulative) distribution function *J>(i) so defined: 

cf>(t) := |$(t) , := P (T < t) = jf* 0(0 dt' . (1.2) 

When the nonnegative random variable represents the lifetime of technical 
systems, it is common to refer to $(t) as to the failure probability and to 

POO 

*(t):=P(T>t)= J (f>(t') dt' = 1 - $(i) , (1.3) 

as to the survival probability, because &(t) and ^(t) are the respective 
probabilities that the system does or does not fail in (0, T]. A relevant 
quantity is the counting function N(t) defined as 

N(t) := max {k\t k < t, k = 0, 1, 2, . . .} , (1.4) 

that represents the effective number of events before or at instant t. In 
particular we have ^f(t) = P(N(t) = 0) . Continuing in the general theory 
we set Fi(t) = 3>(i), fi(t) = (fr(t), and in general 

F k (t) :=P{t k =T l + ...+T k <t) , f k (t) = -F k {t) , k > 1 , (1.5) 
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thus F k (t) represents the probability that the sum of the first k waiting 
times is less or equal t and fk(t) its density. Then, for any fixed k > 1 the 
normalization condition for F k (t) is fulfilled because 

lim F k (t) = P (t k = Ti + . . . + T k < oo) = 1 . (1.6) 

t^oo 

In fact, the sum of k random variables each of which is finite with probability 
1 is finite with probability 1 itself. By setting for consistency Fo(t) = 1 and 
fo(t) = 8(t), the Dirac delta functions, we also note that for k > we have 

P (N(t) =k):=P(t k <t, t k+1 >£)=[* f k (t') (t - t') dt' . (1.7) 

Jo 

We now find it convenient to introduce the simplified * notation for the 
Laplace convolution between two causal well-behaved (generalized) functions 
f(t) and g(t) 

f f(t') g(t - dt' = {f *g) (t) = (g * /) (t) = f /(* - dt' . 
Jo Jo 

Being f k (t) the pdf of the sum of the k iid random variables T\ , . . . , T k with 
pdf 4>(t) , we easily recognize that f k (t) turns out to be the A;- fold convolution 
of 4>{t) with itself, 

fk(t) = {</>*") (t), (1.8) 

so Eq. (1.7) simply reads: 

P (N{t) = k) = ((j)* k * * ) (t) . (1.9) 

Because of the presence of Laplace convolutions a renewal process is suited 
for the Laplace transform method. Throughout this paper we will denote 
by f(s) the Laplace transform of a sufficiently well-behaved (generalized) 
function f(t) according to 

C{f(t);s} = f(s)= e~ st f(t)dt, s>s , 
Jo 

x We find it convenient to recall the formal representation of this generalized function 
in R+ , 

m-.= ^ y *>0. 
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and for 5(t) consistently we will have 5(s) = 1 . Note that for our purposes 
we agree to take s real. We recognize that (1.9) reads in the Laplace domain 



C{P (N(t) = k);s} = U( S )1 k (a) , (1.10) 



where, using (1.3), 



*( a) = IzM. (i.ii) 



2 The Poisson process as a renewal process 

The most celebrated renewal process is the Poisson process characterized by 
a waiting time pdf of exponential type, 

(j)(t) = Ae~ A * , A>0, t>0. (2.1) 

The process has no memory. Its moments turn out to be 

(T) = \, P*> = ^. •■■> ( T ") = ^r, .... (2-2) 

and the survival probability is 

tf(t) :=P{T>t) =e~ xt , t>0. (2.3) 

We know that the probability that k events occur in the interval of length t 
is ^ 

P(N(t) = k) = ^-e- xt , t>0, k = 0,1,2,... . (2.4) 
fc! 

The probability distribution related to the sum of k iid exponential random 
variables is known to be the so-called Erlang distribution (of order k). The 
corresponding density (the Erlang pdf) is thus 

= A 7T^TTi e " At . k = l,2,..., (2.5) 

(ft - 1)! 

so that the Erlang distribution function of order fc turns out to be 
F k (t) = f f k {t') dt' = l- J2 Wf- e" A ' = f) Wf- e-* , t > . (2.6) 

- 70 n=0 n=fc n - 

In the limiting case fc = we recover fo(t) = S(t), Fo(t) = 1, t > 0. 
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The results (2.4)-(2.6) can easily obtained by using the technique of 
the Laplace transform sketched in the previous section noting that for the 
Poisson process we have: 

J(S) = _A_. S(.) = s i 7 , (2.7) 

and for the Erlang distribution: 

f k (s) = m) k = -^, h(s)= [ ^ = -^— k . (2.8) 

(A + S) K S S{A + S) K 

We also recall that the survival probability for the Poisson renewal 
process obeys the ordinary differential equation (of relaxation type) 

^tf(t) = -A*(t), t>0; tf(0+) = l. (2.9) 

3 The renewal process of Mittag-Leffler type 

A "fractional" generalization of the Poisson renewal process is simply 
obtained by generalizing the differential equation (2.9) replacing there the 
first derivative with the integro-differential operator tD* that is interpreted 
as the fractional derivative of order (3 in Caputo's sense, see Appendix. We 
write, taking for simplicity A = 1, 

= -*(t) , * >0, < /3 < 1; *(0+) = l. (3.1) 

We also allow the limiting case (3 = 1 where all the results of the previous 
section (with A = 1) are expected to be recovered. 

For our purpose we need to recall the Mittag-Leffler function as 
the natural "fractional" generalization of the exponential function, that 
characterizes the Poisson process. The Mittag-Leffler function of parameter 
(3 is defined in the complex plane by the power series 

oo n 

It turns out to be an entire function of order (3 which reduces for (3 = 1 to 
exp(z) . For detailed information on the Mittag-Leffler- type functions and 
their Laplace transforms the reader may consult e.g. [4l l8| IT?). 
The solution of Eq. (3.1) is known to be, see e.g. [1, 8, 13J, 

^(t)=Ep(-t f} ), t>0, 0<(3<l, (3.3) 



5 



so 

t(t):=-jV(t) = -j t Ep(-tt 3 ), t>0, 0<(3<1. (3.4) 
Then, the corresponding Laplace transforms read 

^ S ) = rr-F' ^ s ) = r^' °<^<!- ( 3 - 5 ) 

1 + sp 1 + 

Hereafter, we find it convenient to summarize the most relevant features of 
the functions ^(t) and 4>(t) when < (3 < 1 . We begin to quote their series 
expansions for t — > + and asymptotics for i — > oo, 

^ ) = ^ n ( - 1} rcsn + i) ~ — (3 - 6) 

n=0 v ' 

and 

rhm 1 tPU sin(/3vr) + 1) 

In contrast to the Poissonian case (3 = 1, in the case < (3 < 1 for large i 
the functions ^f(t) and <^(i) no longer decay exponentially but algebraically. 
As a consequence of the power-law asymptotics the process turns be no 
longer Markovian but of long-memory type. However, we recognize that 
for < [3 < 1 both functions ^(t), <f>(t) keep the "completely monotonic" 
character of the Poissonian case. Complete monotonicity of the functions 
^f(t) and <p(t) means 

(_l)»_tf(i)>0, (-l) n — 0(i) >0, 71 = 0,1,2,..., t>0, (3.8) 

or equivalently, their representability as real Laplace transforms of non- 
negative generalized functions (or measures), see e.g. [8]. 

For the generalizations of Eqs (2.4) and (2.5)-(2.6), characteristic of the 
Poisson and Erlang distributions respectively, we must point out the Laplace 
transform pair 

C{iP k Ef\-lPy, a }= /3>0, k = 0,1,2,... , (3.9) 

with Ef\z) := —^E^z) , that can be deduced from the book by Podlubny, 
see (1.80) in [T7j. Then, by using the Laplace transform pairs (3.5) and Eqs 
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(3.3), (3.4), (3.9) in Eqs (1.8) and (1.9), we have the generalized Poisson 
distribution, 

P{N{t) = k) = t^Ef{-lP), A; = 0,1,2,... (3.10) 

and the generalized Erlang pdf (of order k > 1), 

fkp-l ... 

The generalized Erlang distribution function turns out to be 
E k (t) = / f k (t') dt' = l-J2— B$\-#) = Y.— E p (- tP ) ■ (3-12) 

JO n U - 7 W H 

n=0 n=fc 

4 The Mittag-Leffler distribution as limit for 
thinned renewal processes 

The procedure of thinning (or rarefaction) for a generic renewal process 
(characterized by a generic random sequence of waiting times {Ifc}) has 
been considered and investigated by Gnedenko and Kovalenko [6]. It means 
that for each positive index k a decision is made: the event is deleted with 
probability p or it is maintained with probability q = 1 — p, with < q < 1. 
For this thinned or rarefied renewal process we shall hereafter revisit and 
complement the results available in [6] . We begin to rescale the time variable 
t replacing it by t/r, with a parameter r on which we will dispose later. 
Denoting, like in (1.5), by F^it) the probability distribution function of the 
sum of k waiting times and by fk(t) its density, we have recursively, in view 
of (1.8), 

/i(t) = 0(t), /*(*)= £ f k -i(t-l/)<t>(t')dt' =(0* fc )(i), k>2. (4.1) 

Let us denote by (Tq jr f)(t) the waiting time density in the thinned and 
rescaled process from one event to the next. Observing that after a 
maintained event the next one of the original process is kept with probability 
q but dropped in favour of the second next with probability pq and, 
generally, n—1 events are dropped in favour of the n-th next with probability 
p' 11 ^ 1 q , we arrive at the formula 

oo 

(T q , r f)(t) = Y,1P n ~ 1 Ut/r)/r . (4.2) 

71=1 
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Let fn(s) = /o°° e ~ st fn(t) dt be the Laplace transform of f n (t). Recalling 
fi(t) = 4>(t) we set fi(s) = cj>(s). Then f n (t/r)/r has the transform f n (rs) = 
(d>(rs)^ , and we obtain (in view of p = 1 — q) the formula 

CT 9 .r?)W = £ IP- 1 (?(«))" = i f q *™ - , (4.3) 
n=l 1 - (1 - ?) 

from which by Laplace inversion we can, in principle, construct the 
transformed process. 

We now imagine stronger and stronger rarefaction (infinite thinning) by 
considering a scale of processes with the parameters r = 5 and q = e tending 
to zero under a scaling relation e = e(5) yet to be specified. Gnedenko and 
Kovalenko have, among other things, shown that if the condition 

0(s) = 1 - a {s) aP + o (a(s) , for s -» 0+ , (4.4) 

where a(s) is a slowly varying function for s — > (H, is satisfied, then we have 
with e = e(5) = a(S) 5^ for every fixed s > the limit relation 

Km 6W ^j =-\, 0<P<1. (4.5) 

5-0 1 _ (1 _ e ) 0(J S ) l + 

This condition is met with a(s) = AM(l/s) if the waiting time T obeys a 
power law with index /3, in the sense of Master Lemma 2 by Gorenflo and 
Abdel-Rehim [7j. The function M(y) is the same as in Master Lemma 2, 
so it varies slowly at infinity, whence M(l/s) varies slowly at zero. The 
proof of (4.5) is by straightforward calculation. Observe the slow variation 
property of a(s) and note that terms small of higher order become negligible 
in the limit. By the continuity theorem for Laplace transforms, see Feller 
[5], we now recognize <po{t) as the limiting density, which we identify, in view 
of (3.2)-(3.5), 

M t) = -± Ef3 (-tP). (4.6) 

So the limiting waiting time density is the so-called Mittag-Leffler density, 
that in the special case (3=1 reduces to the well-known exponential density. 



2 Definition: We call a (measurable) positive function a(y), denned in a right 
neighbourhood of zero, slowly varying at zero if a(cy)/a(y) — > 1 with y — > for every c > 0. 
We call a (measurable) positive function b(y), defined in a neighbourhood of infinity, slowly 
varying at infinity if b(cy)/a(y) — » 1 with y — + oo for every c > 0. Examples: (logy) 7 with 
7 G R and exp (logy/log logy). 
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It should be noted that Gnedenko and Kovalenko in the sixties failed to 
recognize <j)o(s) as Laplace transform of a Mittag-Leffler type function^]. 

5 Renewal processes with reward: the fractional 
master equation and its solution 

The renewal process can be accompanied by reward that means that at 
every renewal instant a space-like variable makes a random jump from its 
previous position to a new point in "space" . "Space" is here meant in a very 
general sense. In the insurance business, e.g., the renewal points are instants 
where the company receives a payment or must give away money to some 
claim of a customer, so space is money. In such process occurring in time 
and in space, also referred to as compound renewal process, the probability 
distribution of jump widths is as relevant as that of waiting times. 

Let us denote by X n the jumps occurring at instants t n , n = 1,2,3, ... . 
Let us assume that X n are iid (real, not necessarily positive) random 
variables with probability density w(x), independent of the waiting time 
density 4>(t). In a physical context the X n s represent the jumps of a diffusing 
particle (the walker), and the resulting random walk model is known as 
continuous time random walk (abbreviated as CTRW) in that the waiting 
time is assumed to be a continuous random variabkfl 

3 Although the Mittag-Leffler function was introduced by the Swedish mathematician 
G. Mittag-Leffler in the first years of the twentieth century, it lived for long time in 
isolation as Cinderella. The term Cinderella function was used in the fifties by the Italian 
mathematician F.G. Tricomi for the incomplete gamma function. In recent years the 
Mittag-Leffler function is gaining more and more popularity in view of the increasing 
applications of the fractional calculus and is classified as 33E12 in the Mathematics Subject 
Classification 2000. 

4 The name CTRW became popular in physics after that in the 1960s Montroll, Weiss 
and Scher (just to cite the pioneers) published a celebrated series of papers on random 
walks to model diffusion processes on lattices, see e.g. [22] and references therein. CTRWs 
are rather good and general phenomenological models for diffusion, including anomalous 
diffusion, provided that the resting time of the walker is much greater than the time it 
takes to make a jump. In fact, in the formalism, jumps are instantaneous. In more recent 
times, CTRWs were applied back to economics and finance by Hilfer [10] , by the authors 
of the present paper with M. Raberto [201 El El HE] > and, later, by Weiss and co-workers 
[15] .However, it should be noted that the idea of combining a stochastic process for waiting 
times between two consecutive events and another stochastic process which associates a 
reward or a claim to each event dates back at least to the first half of the twentieth century 
with the so-called Cramer-Lundberg model for insurance risk, see for a review [5]. In a 
probabilistic framework, we now find it more appropriate to refer to all these processes as 
to compound renewal processes. 
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The position x of the walker at time t is 



x(t) 



N(t) 

;(0) + £X fc . 



(5.1) 



k=l 



Let us now denote by p(x, t) the probability density of finding the random 
walker at the position x at time instant t . We assume the initial condition 
p(x,0) = S(x) , meaning that the walker is initially at the origin, x(0) = . 
We look for the evolution equation for p(x, t) of the compound renewal 
process. Based upon the previous probabilistic arguments we arrive at 



p(x,t) =6(x)V(t)+ / <t>{t -t') 



w(x — x) p{x , t) dx 



dtl , (5.2) 



called the integral equation of the CTRW. In fact, from Eq. (5.2) we 
recognize the role of the survival probability vE'(t) and of the densities 
4>{t) , w(x) . The first term in the RHS of (5.2) expresses the persistence 
(whose strength decreases with increasing time) of the initial position x = 0. 
The second term (a space-time convolution) gives the contribution to p(x, t) 
from the walker sitting in point x' £ R at instant t' < t jumping to point x 
just at instant t , after stopping (or waiting) time t — t' . 

The integral equation (5.2) can be solved by using the machinery of the 
transforms of Laplace and Fourier. Having introduced the notation for the 
Laplace transform in sec. 1, we now quote our notation for the Fourier 
transform, JF{/(x);k} = f(n) = e lKX f(x) dx (k E R), and for the 

corresponding Fourier convolution between (generalized) functions 



(/i * h) (x) 



+oo 



h{x')f 2 {x-x')dx'. 



Then, applying the transforms of Fourier and Laplace in succession to 
(5.2) and using the well-known operational rules, we arrive at the famous 
Montroll- Weiss equation, see [16] . 



p(k,s) 



1 — 4>(s) w(k) 



(5.3) 



As pointed out in [7], this equation can alternatively be derived from 
the Cox formula, see [2] chapter 8 formula (4), describing the process as 
subordination of a random walk to a renewal process. 
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By inverting the transforms one can, in principle, find the evolution 
p(x,t) of the sojourn density for time t running from zero to infinity. In 
fact, recalling that |u5(k)| < 1 and \<j>{s)\ < 1, if k ^ and s / 0, Eq. (5.3) 
becomes 

oo 

K« jS ) = *( s ) £[0(*)^)] fc ; (5-4) 

k=0 

this gives, inverting the Fourier and the Laplace transforms and taking into 
account Eqs. (1.9)-(1.10), 



p(x 



t) = J2P(N(t) = k)w k {x), (5.5) 



k=0 



where Wk(x) = (w* k ^J (x), in particular wo(x) = 5(x), w\{x) = w(x). 

A special case of the integral equation (5.2) is obtained for the compound 
Poisson process where (f)(t) = e _i (as in (2.1) with A = 1 for simplicity). 
Then, the corresponding equation reduces after some manipulations, that 
best are carried out in the Laplace-Fourier domain, to the Kolmogorov- Feller 
equation: 

d r +ao 

— p(x,t) = —p(x,t) + / w{x — x')p(x', t) dx' , (5-6) 

Ot J —oo 

which is the master equation of the compound Poisson process. In this case, 
in view of Eqs (2.4) and (5.5) the solution reads 

oo 

p{x,t) = Y J Tr,e~ t Wk{x). (5.7) 

k=Q K - 

When the survival probability is the Mittag-Lemer function introduced in 
(3.3), the master equation for the corresponding fractional version of the 
compound process can be shown to be 

/+oo 
w(x - x')p(x',t)dx' , Q</3<1, (5.8) 
-oo 

where tD* denotes the time fractional derivative of order [3 in the Caputo 
sense. For a (detailed) derivation of Eq (5.8 ) we refer to the paper by 
Mainardi et al. [2], in which the results have been obtained by an approach 
independent from that adopted in a previous paper by Hilfer and Anton [TTj . 
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In this case, in view of Eqs (3.10) and (5.5), the solution of the fractional 
master equation (5.8) reads: 



In [9] we have, under a power law regime for the jumps, investigated for 
Eq. (5.8) the so-called diffusive or hydrodynamic limit, obtained by making 
smaller all jumps by a positive factor h and accelerating the process by 
a large factor properly related to h, then letting h tend to zero. In this 
limit the master equation (5.8) reduces to a space-time fractional diffusion 
equation. This is also the topic of the recent paper by Scalas et al. [21] and, 
in a more general framework, of the paper by Gorenflo and Abdel-Rehim 



Conclusions 

We have provided a fractional generalization of the Poisson renewal processes 
by replacing the first time derivative in the relaxation equation of the 
survival probability by a fractional derivative of order (0 < /3 < 1). 
Consequently, we have obtained for < (3 < 1 non-Markovian renewal 
processes where, essentially, the exponential probability densities, typical 
for the Poisson processes, are replaced by functions of Mittag-Leffler type, 
that decay in a power law manner with an exponent related to (5. 

The distributions obtained by considering the sum of k iid random vari- 
ables distributed according to the Mittag-Leffler law provide the fractional 
generalization of the corresponding Erlang distributions. Furthermore, 
the Mittag-Leffler probability distribution is shown to be the limiting 
distribution for the thinning procedure of a generic renewal process with 
waiting time density of power law character. 

Then, our theory has been applied to renewal processes with reward, so 
can be considered as the fractional generalization of the compound Poisson 
processes. In such processes, occurring in time and in space, also the 
probability distribution of the jump widths is relevant. The stochastic 
evolution of the space variable in time is modelled by an integro-differential 
equation (the master equation) which, by containing a time fractional 
derivative, can be considered as the fractional generalization of the classical 
Kolmogorov-Feller equation of the compound Poisson process. For this 
master equation we have provided the analytical solution in terms of iterated 
derivatives of a Mittag-Leffler function. 




k=0 



(5.9) 
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Appendix: The Caputo fractional derivative 



The Caputo fractional derivative provides a fractional generalization of the 
first derivative through the following rule in the Laplace transform domain, 



c{ t D^f(t);s} = s' 3 f(s)-s' 3 - 1 f(0 + ), O<0<1, s>0. 
hence turns out to be defined as, see e.g. [HE], 



tD?f(t) 



/ (1) (r) 

T(i-p)j Q (t- T y 



dr . 



0<P<1, 
= 1. 



It can alternatively be written in the form 



tD?f(t) 
1 



d 



fir) 



dr 



r(l-/3) dtJo (t-r)/» r(i-/3) 



/(o^ 



d /•t/(r)-/(0+) 



r i 



dr , < /3 < 1 



(Al) 



(A2) 



(A3) 



0) dtVo (t-r)/ 3 

The Caputo derivative has been indexed with * in order to distinguish it 
from the classical Riemann-Liouville fractional derivative t-D^, the first term 
at the R.H.S. of the first equality in (A. 3). As it can be noted from the last 
equality in (A. 3), the Caputo derivative provides a sort of regularization at 
t = of the Riemann-Liouville derivative; for more details see [8]. 
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